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We present an ab initio calculation of the shear viscosity as a function of interaction strength in a 
two-component unpolarized Fermi gas near the unitary limit, within a finite temperature quantum 
Monte Carlo (QMC) framework and using the Kubo linear-response formalism. The shear viscosity 
decreases as we tune the interaction strength from the Bardeen-Cooper-Schrieffer side of the 
Feshbach resonance towards Bose-Einstein condensation limit and it acquires the smallest value for 
—^ ~ 0.4, with a minimum value of - « 0.2 jS-, which is about twice as small as the value 

reported for experiments in quark-gluon plasma f Ig^p ^ The Fermi gas near unitarity thus 

emerges as the most “perfect fluid” observed so far in nature. The clouds of dilute Fermi gas near 
unitarity exhibit the unusual attribute that, for the sizes realized so far in the laboratory or larger 
(less than 10® atoms), can sustain quantum turbulence below the critical temperature, but at the 
same time the classical turbulence is suppressed in the normal phase. 

PACS numbers: 03.75.Ss, 05.30.Fk, 05.60.Gg, 51.20.-l-d 


I. INTRODUCTION 


The ultracold atoms provide an ideal laboratory for 
very precise experimental and theoretical studies of an 
enormous range of quantum mechanical phenomena. A 
large set of studies focused on interacting Fermi gases at 
unitarity, where such systems exhibit remarkable prop¬ 
erties. The average inter-particle separation Ri 
is large compared to the effective range of interaction 
Tgjf, but small compared to the scattering length |a|, i.e., 
0 T- <C 1 <C -T oo (where n is particle 

density). In this limit these systems acquire universal 
properties and they have been widely studied over the 
last dozen years both experimentally and theoretically 
(see reviews [11-13 )■ ^ series of recent experiments (B-jll 
revealed the nearly ideal hydrodynamic behavior of the 
resonantly interacting Fermi gas, characterized by a very 
low shear viscosity coefficient to the entropy density ra¬ 
tio - very close to the conjectured bound originating from 
holographic duality methods @,[13 , called KSS bound. 
An “ideal fluid” which follows the laws of ideal hydrody¬ 
namic^ is a fluid in which dissipative processes are ab¬ 
sent ll|. The superfluid component of a quantum fluid 


below the critical temperature is treated as a physical 
realization of an ideal fluid [13 ■ The normal component 
is characterized by a finite viscosity, which in an infinite 
medium tends to infinity at absolute zero temperature, 
due to the contributions of phonons [l3- In a finite sys¬ 
tem however the longest phonon wavelength is of the or¬ 
der of the size of the container and with the number of 
excited phonons oc this has the consequence that the 
viscosity, while approaching the zero temperature, never 
formally diverges and a liquid is effectively in a collision¬ 
less regime. 

A hypothetical physical system which saturates the 


KSS bound and has the lowest possible value of the shear 
viscosity is often referred to as “the perfect fluid”. A re¬ 
cent measurement of the shear viscosity of a unitary gas 
close to unitarity [3 provides the value ^ « 0.5^, while 
KSS bound is ~ 0.08^. Thus the unitary Fermi 

gas (UFG) appears as the system being very close to 
“perfectness.” Another system is the quark-gluon plasma 
created in relativistic heavy ion collisions 11311 , where the 
value f Iqqp has been reported [lj|- 

Simultaneously, an impressive effort has been made in 
order to theoretically determine transport coefficients of 
the UFG. A plethora of theoretical methods has been 
used to estimate these transport coefficients, both for 
homogeneous and trapped systems [isl - l^ . Typically 
these theoretical predictions differ both quantitatively 
and qualitatively. Among the methods used, only the 
works [13 Hi present ab initio calculations, obtained 
within the very powerful quantum Monte Garlo (QMG) 
framework, where the errors can be quantified quite ac¬ 
curately. In these papers, however, the shear viscosity 
has been evaluated only at unitarity. In Fig. [T]we com¬ 
pare these calculations with recent experimental data 
for a uniform system extracted from measurements for 
trapped systems [1|- Agreement for absolute value of 
shear viscosity has not been obtained, which however 
does not rule out the QMC results. However, we observe 
a very good qualitative agreement. Experimentally, the 
results for uniform systems @] have been inferred from 
data for trapped systems under a number of assumptions 
(see Supplement of Ref. @|)- To what extent these as¬ 
sumptions are valid, especially in the most interesting 
low temperature regime, is not clear. In particular, vis¬ 
cous hydrodynamics has been used in Ref. @| to describe 
the full dynamics of the cloud, even at the periphery of 
the cloud, where the density is very low and the colli- 
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FIG. 1: (Color online) Comparison of the QMC predictions 
taken from Ref. [IJ (blue points) with the results extracted 
from experimental data for trapped system (red points) [^. 
Horizontal and vertical axes show the temperature in units of 

Fermi energy ep = 2 m dimensionless shear viscos¬ 

ity rj/n, respectively. The theoretical results are provided for 
different lattice sizes = 8,10 and 12. The scatter of these 
points is a measure of the numerical accuracy of the available 
QMC predictions. The solid black line shows kinetic theory 
/ ™ \ 3/2 

prediction ^ = 2.77 ( ^ ) ■ The phonon contribution to 

the viscosity, which is not accounted for in QMC calculations 
in a finite volume and evaluated in Ref. El , is shown as 
a dot-dashed (brown) line. The vertical black dotted lines in¬ 
dicate the critical temperature of superfluid-to-normal phase 
transition for Tc = 0.15(1) and the onset of Cooper pair 
formation, T* ~ 0.22 £f, resp ectively - both extracted from 
QMC simulations in Refs and [2^. 


sion rate is strongly suppressed and the system is in the 
collisionless regime [2l. In both cases we observe the 
rapid decrease of the shear viscosity well above critical 
temperature of superfluid-normal phase transition Tc- It 
is notable that this feature is present only in the QMC 
and the pseudogap T-matrix theory El predictions. 

Recently it has been reported that the system is even 
closer to the prefect fluid limit if it is slightly beyond 
the unitary point . In this experiment it was revealed 
that the lowest viscosity coefficient can be obtained for an 
interaction strength corresponding to ~ 0.25, where 
kp = is Fermi momentum. Since the QMC 

method provides results which are at least in very good 
qualitative agreement with experiment for the unitary 
limit (i.e., —>■ 0), in this work we extend our studies 

beyond unitarity point. 


II. METHOD 


In order to determine the shear viscosity coefficient of 
the ultracold atomic gas we employ the QMC technique 
with auxiliary fields on the lattice, which provides numer¬ 


ical results with controllable accuracy, up to quantifiable 
systematic uncertainties (for details see Ref. [27| ). These 
simulations are very similar to those performed at unitar¬ 
ity when |a| = 00 (Refs. jUjUS])- Here we briefly describe 
the main aspects of the computational process, focusing 
mainly on required modifications in order to study the 
system beyond the unitarity. Henceforth we use the sys¬ 
tem of units: h = m = kp = 

For simulations we employed a cubic lattice of size 
Nx = Ny = Nz = 10, with lattice spacing I = 1 and 
average number densities n ~ 0.04. As shown in the 
Supplemental Material of Refs. [1^ [2^ the systematic 
errors are no more than 10% for this lattice size and the 
errors are related mainly to corrections coming from the 
nonzero effective range and from the exclusion of the 
universal high momenta tail in the occupation probabil¬ 
ity due to the finite momentum cut-off — J- 

the context of computation of transport coefficient it is 
important to have as small as possible statistical errors, 
so as to minimize the errors arising from the analytical 
continuation from the imaginary time to real frequencies 
as explained below. Thus we generate an ensemble con¬ 
taining about 10"^ uncorrelated samples in order to get 
statistical accuracy below 1%. 

In order to extract the shear viscosity within a QMC 
framework one calculates the imaginary-time (Euclidean) 
stress tensor-stress tensor correlator 

Gn(g,r) = :i(fl^-^)(T)nLf (0)), (I) 

at zero momentum q = 0. As shown in Ref. [isj ]. for zero- 
range interaction it is sufficient to use only the kinetic 
part of the stress tensor: 

P^Pv^xiP)^x(P)- (2) 

pA=t,4- 


The average is performed in the grand canonical ensem¬ 
ble, at fixed temperature T = ^ and chemical poten¬ 
tial, = exp[T(7? — exp[— t(IT — fiN)], 

where H is the Hamiltonian of the system, /r is the chem¬ 
ical potential, and N is the particle number operator. In 
order to capture the physics of a dilute fermionic gas 
it is sufficient to use a zero-range two-body interaction 
V{ri — r 2 ) = —gS{ri — r 2 ), where the coupling constant 
g can be tuned to fix the value of the s-wave scatter¬ 
ing length a using a standard renormalization pro¬ 
cedure of the coupling constant on the lattice. After this 
procedure the contact interaction acquires a finite effec¬ 
tive range - for the parameters of presented simulation 
kpr^g ~ 0.43. In this case the interaction part of the 
stress tensor also contributes to the correlator. How¬ 
ever, we relegate all the corrections arising from the fi¬ 
nite effective range to the systematic errors, which were 
estimated to be less than 10% for used lattice (see the 
Supplemental Material of Ref. 2^ for extensive discus¬ 
sion in the context of lattice computation of transport 
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coefficients). It is legitimate to drop this contribution, 
particularly since in the limit kpr^g —>■ 0 it is vanish¬ 
ing anyway, a fact also consistent with previous QMC 
studies performed at unitarity for different densities. We 
emphasize that the dominant source of uncertainties is 
introduced by the analytic continuation procedure (de¬ 
scribed next). These are of the order 10%, as one can 
judge from Fig. [T 1 where the error bars show contribu¬ 
tion to the uncertainties only from this source. 

The frequency dependent shear viscosity 77 ( 0 ;) is ex¬ 
tracted via the analytic continuation of the imaginary- 
time correlator to real frequencies. This procedure is 
equivalent to solving the integral equation: 


Gii{q = 0 ,r) = - 

TT 


’ cosh[a;(T-^/2)] 

sinh(c.g/2) 


The static shear viscosity rj is defined as ry = limj^^o 
The correlator Gn is sampled only for a finite set of points 
and in a finite imaginary-time interval, and its evaluation 
is affected by the statistical noise, which we minimize by 
using a quite high statistics. This integral equation ([3]) 
belongs to a class of numerically ill-posed problems. 
Therefore, the use of special techniques is warranted in 
order to extract numerically stable results. We have em¬ 
ployed an approach which combines two complementary 
methods: the Singular Value Decomposition (SVD) and 
the self-consistent Maximum Entropy Method (MEM), 
both described in great detail in Ref. [1^ and in the 
Supplemental Material of Refs. [ 2 ^ 281. The stabiliza¬ 


tion procedure requires a priori information about the 
solution 77 ( 0 ;). The a priori information used is re¬ 
lated to the known properties: the non-negativity of the 
shear viscosity 77 ( 0 ;) ^ 0 , the asymptotic tail behavior 
77 ( 0 ; —>• 00 ) = 0 


subsequent corrections [l 8 l l31jl 


Ifd.. 

TT Jo 


77 ( 0 ;) 



£ C 
3 ~ 12™’ 


(4) 


where C is Tan’s contact density [s^ , and £ is the energy 
density. Both quantities are obtained within the same 
QMC simulation. The contact density C was extracted 
from the analysis of the tail of the numerically gener¬ 
ated momentum distribution, which for sufficiently large 
momenta decays as n{p) ^ a similar technique was 
used in Ref. [s^ and subsequent studies. Moreover, self- 
consistent MEM requires an appropriately chosen class 
of a priori models for the solution. Based on our past 
experience for the unitary limit (^ . , we determined 

that the expected suitable models for the frequency de¬ 
pendent shear viscosity 77 ( 0 ;) consist of Lorentzian-like 
structures at low frequencies, smoothly evolving into the 
asymptotic tail behavior: 


+ [l-/(w: {“n « 2 })] c}), (5) 


where 


f{uj, {ofi, a 2 }) = exp(-aia 2 ) 


exp(a]^a;) — 1 
1 -I- exp(ai(w — a^)) 


( 6 ) 


and 


J^(w,{m,7,c}) = C-7--2- 

TT (w — -I- y"’ 

The parameters {/i, 7 , c, 0 : 2 } describe admissible de¬ 
grees of freedom of the model and are adjusted automat¬ 
ically in a self-consistent manner. 


III. RESULTS 


In Fig. [21 the dimensionless static shear viscosity ^ 
is shown for three selected temperatures ^ = 0.26, 0.4 
and 0.5 as a function of interaction strength \/{akp). 
Only one point, for ^ = 0.26 and = 0.2, cor¬ 
responds to the system being in the superfluid phase, 
while all other points correspond to the system either 
in the normal phase or in the “pseudogap” regime. We 
emphasize that our QMC simulations are fully consis¬ 
tent with the existence of the “pseudogap” regime, i.e., 
a temperature regime above the critical temperature Tc 
where many Cooper pairs are present, even though the 
superfluidity is lost 2^, 34 1 • The presence of pairs above 
the critical temperature is a property of the UFG well 
established in ab initio calculations. Naturally, this is 
also beyond controversy in the BEC limit as well, where 
the critical temperature is well below the dimer binding 
energy, Tc The presence of a pseudogap regime 

is to some extent also confirmed by experiments |35l - l37j |. 
The present QMC results clearly show that the static 
shear viscosity decreases as we enter into BEC regime. 
As we tune the system towards the BEC side of the Fes- 
hbach resonance the number of Monte Carlo samples re¬ 
quired to get an acceptable signal-to-noise ratio increases 
rapidly. For this reason we were not able to perform re¬ 
liable inversions of the Eq. ([3|) for interaction strengths 
beyond > 0.2 for temperature ^ = 0.26 and be¬ 
yond > 0.5 for temperatures ^ = 0.4 and 0.5. The 
ratio of the shear viscosity to the number density — ap¬ 
pears to have a relatively weak temperature dependence 
above Tc- This behavior appears to be confirmed by the 
experiment [2^ as well. Our results for temperatures 
^ = 0.4 and 0.5 suggest that there exists a minimum 
for the static shear viscosity located on the BEC side of 
resonance for an interaction strength corresponding to 
Ri 0.4. However, the present QMC uncertainties do 
not permit an accurate determination of the minimum 
position. 

Comparing our results with experimental data [ 2 ^ we 
find qualitative agreement. While the experiment re¬ 
ports existence of minimum for the shear viscosity for 
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FIG. 2: (Color online) The dimensionless static shear vis¬ 
cosity ^ as a function of interaction strength for a 10^ 
lattice for selected temperatures: T = 0.26eir - solid (red) cir¬ 
cles, T = 0.4e_F - (blue) squares and T = O.Sef - (green) dia¬ 
monds. The error bars show contribution to the computation 
uncertainty generated by the analytic continuation procedure 
only. 


~ 0.25, the ab initio prediction sets it for slightly 
higher strengths. In general, the minimum position can 
be temperature dependent, but experimental data sug¬ 
gest that this dependence is rather weak, in apparent 
agreement with our findings. Note, however, that the 
measurements are performed in trapped systems and only 
the trap-averaged viscosity (^) = f r]{r)d^r is re¬ 
ported in Ref. [26|. The reduced temperature ^ as 
well as the interaction strength are position depen¬ 
dent in a trap and they diverge to infinity as one ap¬ 
proaches the trap edges, where the gas enters the colli¬ 
sionless regime and where the hydrodynamical approach 
is inapplicable The trapped experiment probes all 
reduced temperatures and interaction strengths, starting 
from the values at the center of the cloud up to very 
high values. Moreover, the regions further from the cen¬ 
ter of the trap also contribute with increasingly higher 
weights in such averaged quantities, masking to a large 
extent the information about the inner regions of the trap 
and making its evaluation very challenging. To what ex¬ 
tent the averaging procedure affects the results for shear 
viscosity is not clear at this time, and since a reliable 
validation of the theoretical predictions against the ex¬ 
perimental data requires knowledge of the shear viscos¬ 
ity for all temperatures and interaction strengths, such a 
comparison is beyond the scope of the present work. A 
similar disagreement with experimental values was noted 
in the analysis performed by Bluhm and Schafer . For 
example, at unitarity the values extracted in the experi¬ 
mental analysis of Ref. Q of the ratio ^ exceed unity at 
temperatures ^ > 0.3. The calculated ratio ^ within 
the kinetic theory [s^ attains such values only for signif¬ 


icantly larger temperatures ^ Ri 1 at unitarity (see Fig. 
[1} and even greater temperatures for positive values of 
the scattering length. The kinetic theory is appropriate 
for temperatures above T*, where pairs have completely 
dissociated and where the collision integral is more or less 
well defined. Below T* a more complex kinetic approach 
is required, which should include dimer-dimer, fermion- 
fermion, and dimer-fermion collisions as well as a dimer- 
to-two-fermions and its time-reverse processes. On the 
other hand, when extrapolated, the kinetic theory and 
the present QMC results appear, surprisingly to some 
extent, to be in agreement. 

In order to confront the QMC results with the KSS 
conjuncture one has to have information about the en¬ 
tropy density s = y. This can be extracted from static 
observables including the energy E, the chemical po¬ 
tential /i, and the contact C, which are easily obtained 
within the QMC framework. Combining the basic ther¬ 
modynamic relation (where T, P, and N are, respec¬ 
tively, temperature, pressure, and particle number) 


E = TS - PV + ^lN 
together with Tan’s pressure relation (39| 

3F 127raF’ 


one can show that 

S{x,y) _ ^ix,y) - Cix,y) + ^C{x,y)y 
N X 


( 8 ) 

(9) 


( 10 ) 


where we introduced the following dimensionless quan¬ 
tities: the reduced temperature a: = ^, the strength 
of the interaction y = {akF)~^, the Bertsch parameter 
^ reduced chemical potential C = and 

the reduced contact parameter C = 

In Fig. [3] we show the ratio of the shear viscosity to 
the entropy density These results suggest that there 
is minimum for this ratio located at ~ 0.4. How- 
ever, the present QMC uncertainties do not permit one to 
make a very precise determination of its location. More¬ 
over, the QMC data also reveal a rather weak tempera¬ 
ture dependence, similar to experimental findings. 

The - ratio at the minimum -I . ~ 0.2, is about 

2.5 times smaller than its value at the unitary limit, and 
only 2.5 times larger than the KSS holographic bound 
In the inset we also provide the entropy dependence on 
the coupling constant for our selected temperatures. For 
the temperature ^ = 0.26 the value of the entropy for 
= 0.2, which visually appears to deviate from the 
smooth pattern, reflects the fact that there is a phase 
transition from the normal to the superfluid state as we 
increase interaction strength at a fixed temperature. 
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FIG. 3: (Color online) The ratio of the shear viscosity to en¬ 
tropy density ^ as a function of interaction strength for 
a 10® lattice. The notation for theoretical series is identical 
to Fig. (2] The black dashed line indicates KSS bound: 

The error bars encapsulate uncertainty originating from the 
analytic continuation procedure and statistical uncertainty of 
the entropy density determination. The inset shows corre¬ 
sponding values of entropy per particle. 


IV. QUANTUM & CLASSICAL TURBULENCE 


The viscosity is at the root of classical turbulence. In 
an ideal fluid turbulence does not exist, while it can de¬ 
velop in either a normal fluid or in a superfluid. In super¬ 
fluids viscosity was expected to play a minor role (only 
in the normal component) and there was no clear route 
to turbulence when the temperature tends to zero and 
the viscosity therefore decreases. However, in 1955 Feyn¬ 
man [ 4 ^ conjectured that in superfluids the crossing and 
reconnection of quantized vortices could lead to “quan¬ 
tum turbulence,” a field which since then is one of the 
most active areas of research in the physics of liquid he¬ 
lium 3 and 4 Feyman’s conjecture was demon¬ 

strated to be the correct theoretical mechanism in a di¬ 


lute Bose superfluid [44| and in a dilute Fermi superfluid 
as well 45, ka. The theoretical framework of classical 


turbulence laid by Kolomogorov appears to explain 
many features (though not all) of quantum turbulence in 
liquid helium [42j| . At the same time many characteris¬ 
tics, such as the non-Gaussian velocity distributions, are 
drastically different ( 4 ^ in classical and quantum tur¬ 
bulence. The non-Gaussian velocity distributions have 
been predicted to emerge also in the UFG [i^. The 
quark-gluon plasma and the Fermi gas in the unitary 
regime above the critical temperature are two physical 
systems, which are not superfluid, but in which the the 
shear viscosity attains extremely low values. Quantized 
vortices do not exist either in the UFG above Tc or in 
the quark-gluon plasma and the dynamics of these sys¬ 
tems can be very close to that of a hypothetical classical 


fluid with zero viscosity, called the ideal fluid, if exis¬ 
tence of classical turbulence is prohibited. Indeed, this 
is a case of ultracold fermionic gases produced experi¬ 
mentally. Glassical fluid hydrodynamics is governed by 
dimensionless numbers, where the most important is the 
Reynolds number Re = , v and L are characteristic 

velocity and linear dimension describing flow. Classical 
turbulence in three-dimensional systems is achieved for 
values of the Reynolds number of the order of 10'^. 


In a UFG quantized vortices, and therefore quantum 
turbulence, can exist in clouds with as little as 500-1000 
fermions 4^, 4^ and for flow velocities v ~ 0.7uf, which 
are larger than the Landau’s critical velocity Vc ~ QAvp- 
The largest cold atomic clouds created so far in the lab¬ 
oratory have 0(10®) atoms. 


One can estimate the Reynolds number for a UFG us¬ 
ing as characteristic scales: number of atoms V in a 
cloud that defines density n = N/L^, the critical veloc¬ 
ity Vc = QAvf = 0.4 ^*-®^^y — 1.2N^^^h/mL, and the 
minimal value of the shear viscosity ^ « 0.2. With 
this one obtains for the Reynolds number the values 
Re < Ki 60 for N = 1000 and Re Ri 620 for 

iV = 10 ®. One can argue that one can attain higher 
values of the Reynolds number by increasing the flow ve¬ 
locity by a factor 10 so as to reach Re ~ 10“^. In a 
dilute Fermi gas near the unitary point the scattering 
cross section is on average a ~ A^Trjkl^. If the flow veloc¬ 
ity is increased to u ~ IOuf, the cross section decreases 
by a factor of k'^/kp 100 and the mean free path 1/na 
becomes very large, comparable or exceeding the size of 
any atomic cloud created so far in the laboratory and 
the system enters the collisionless regime. One could al¬ 
ternatively contemplate an increase in the linear size of a 
cloud by a factor of 10, thus up to cloud particle numbers 
0(10®) (a size likely difficult to achieve for condensates), 
in order to increase the Reynolds number by an order of 
magnitude. 


The Fermi gas in the unitary regime is thus a rather 
unique physical system; below the critical temperature 
the system is superfluid and can sustain quantum tur¬ 
bulence in rather small clouds, while above the criti¬ 
cal temperature the turbulent dynamics is strongly sup¬ 
pressed for any current experimental realizations for a 
very wide range of flow velocities and cloud sizes. The 
almost “death” of classical turbulence above Tc and its 
revival into a new “body,” the quantum turbulence be¬ 
low Tc, makes the unitary Fermi gas the unitary regime 
a quite unique physical system. The small value of the 
shear viscosity in a dilute Fermi gas near unitarity, which 
is attained in the normal phase, requires that the onset of 
classical turbulence be achieved in relatively large clouds, 
so far not realized experimentally. 


In summary, we have presented ah initio results for 
interaction strength dependence of the static shear vis¬ 
cosity and the shear viscosity to the entropy density ratio. 



















6 


Both quantities decrease as we tune interaction strength 
from the BCS side of the unitarity point towards the BEC 
limit. The results suggest that the Fermi gas in the uni¬ 
tary regime is the closest physical system known to being 
a “perfect fluid,” for an interaction strength correspond¬ 
ing to ~ 0.4, with a minimum value of shear vis¬ 


cosity to entropy density 


~ 0 . 27 ^ about twice as 

in KB 

small than the value reported for the quark-gluon plasma 
f Iqgp ^ simulations qualitatively con¬ 

firm the experimental observation of Ref. [ 2 ^ , that shear 
viscosity attains a minimum on the BEC side of the uni¬ 
tary point, albeit for a stronger value of the coupling 
constant. 
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